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A  fast  parallel  iterative  method  is  proposed  for  the  solution  of 
linear  systems  arizing  from  Finite  Element  discretization  of  the  time 
harmonic  acoustics  of  coupled  fluid-solid  systems  in  fluid  pressure  and  solid 
displacement  formulation.  The  method  generalizes  the  FETI-H  method 
for  the  Helmholtz  equation  to  elastic  scattering.  The  fluid  and  the  solid 
domains  are  decomposed  into  non-overlapping  subdomains.  Continuity  of 
the  solution  enforced  by  Lagrange  multipliers.  The  system  is  augmented 
by  duplicating  the  degrees  of  freedom  on  the  wet  interface.  The  original 
degrees  of  freedom  are  then  eliminated  and  the  resulting  system  is  solved 
by  the  GCR  method  preconditioned  by  a  subspace  correction.  In  each 
iteration,  the  method  requires  the  solution  of  one  independent  acoustic 
problem  per  subdomain,  and  the  solution  of  a  coarse  problem  with  several 
degrees  of  freedom  per  subdomain.  Computational  results  show  that  the 
method  is  scalable  with  the  problem  size,  frequency,  and  the  number  of 
subdomains.  The  number  of  iterations  was  mostly  about  same  as  the 
number  of  iterations  of  the  FETI-H  method  for  the  related  Helmholtz 
problem  with  Neumann  boundary  condition  instead  of  elastic  scatterer. 
Convergence  is  explained  from  the  spectrum  of  the  iteration  operator. 
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1.  INTRODUCTION 

The  method  of  Finite  Element  Tearing  and  Interconnecting  (FETI)  was  proposed 
by  Farhat  and  Roux  [11]  for  solving  linearized  elasticity  problems.  The  FETI 
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method  consists  of  decomposing  the  domain  into  non-overlapping  subdomains. 
Lagrange  multipliers  are  used  to  enforce  that  the  values  of  the  degrees  of  freedom 
coincide  on  the  interfaces  between  the  subdomains  interfaces.  All  original  degrees 
of  freedom  are  then  eliminated,  leaving  a  dual  system  for  the  Lagrange  multipliers, 
which  is  solved  by  preconditioned  conjugate  gradients.  For  details,  further 
developments,  and  theoretical  analysis  of  the  FETI  method  for  positive  definite 
problems,  see  [1,  4,  5,  9,  10,  13,  14,  16,  17,  18,  20,  21]  and  references  therein. 

Variants  of  the  FETI  method  for  the  Helmholtz  equation  of  scattering,  were 
proposed  by  De  La  Bourdonnaye  et  al.  [3],  Farhat  et  al.  [6,  7,  8],  and  further 
developed  by  Tezaur  et  al.  [22].  In  this  paper,  we  present  an  extension  of  the 
method  of  [7,  8],  called  FETI-H,  to  the  case  of  an  elastic  scatterer. 

The  paper  is  organized  as  follows.  In  Section  2,  we  state  the  boundary 
value  problem  of  coupled  fluid-solid  acoustics.  Section  3  reviews  the  variational 
formulation  and  Finite  Element  discretization.  The  discrete  equations  are  put  in 
a  substructured  form  in  Section  4,  which  leads  to  the  iterative  method,  presented  in 
Section  5.  Section  6  contains  computational  results.  An  explanation  of  some  of  the 
convergence  properties  of  the  method  from  the  spectrum  of  the  iteration  operator 
is  given  in  Section  7.  Section  8  is  the  conclusion. 

2.  THE  SCATTERING  PROBLEM 

We  consider  an  acoustic  scattering  problem  with  an  elastic  scatterer  completely 
immersed  in  a  fluid.  Let  fi  and  f2e  be  bounded  domains  in  9?",  =  2,  3,  fie  C  f2, 
and  let  12/  =  12  \  f2e.  Let  v  denote  the  exterior  normal.  Let  <9f2  be  decomposed 
disjoint  subsets,  <9f2  =  U  r„  U  Ta.  The  domain  f2/  is  filled  with  a  fluid.  The 
acoustic  pressure  at  time  t  is  assumed  to  be  of  the  form  R epelult ,  where  p  is  complex 
amplitude  independent  of  t.  The  amplitude  p  is  governed  by  the  Helmholtz  equation 

Ap+k2p  =  0  in  12/ ,  (1) 

with  the  boundary  conditions 

dp  dp 

P  =  P0  onrd,  —  =  0  on  Tn,  —  +  ikp  =  0  on  ro,  (2) 

av  av 

where  k  =  w/c/  is  the  wave  number  and  c/  is  the  speed  of  sound  in  the  fluid. 
The  boundary  conditions  (2)  model  excitation,  sound  hard  boundary,  and  outgoing 
boundary,  respectively.  The  amplitude  of  the  displacement  u  of  the  elastic  body 
occupying  the  domain  f2e  satisfies  the  elastodynamic  equation 

V  •  t  +  u>2  peu  =  0  in  f2e,  (3) 

where  r  is  the  stress  tensor  and  pe  is  the  density  of  the  solid.  For  simplicity,  we 
consider  an  isotropic  homogeneous  material  with 

t  =  XI(y  -u)  +  2pe(u),  ey(u)  = -(— +  ^-),  (4) 

where  A  and  p  are  the  Lame  coefficients  of  the  solid. 
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Let  r  =  dQe  be  the  wet  interface.  On  T,  the  fluid  pressure  and  the  solid 
displacement  satisfy 


u-u  = - r— ,  v  •  t  ■  v  =  —  p,  v  x  t  ■  v  =  0,  (5) 

pfU)Z  ov 

where  pf  is  the  fluid  density.  The  interface  conditions  (5)  model  the  continuity 
of  normal  displacement,  the  balance  of  normal  forces,  and  zero  tangential  tension, 
respectively,  cf.,  e.g.,  [23]. 

3.  VARIATIONAL  FORM  AND  DISCRETIZATION 

We  use  the  following  standard  variational  form  and  discretization  by  conforming 
elements,  cf.,  e.g.,  [19].  Define  the  spaces  Vf  —  {q  £  R1(f If)  \  q  =  0  on  T^}, 
Ve  =  (H1( Oe))n,  where  H 1  is  the  Sobolev  space  of  generalized  functions  with 
square  integrable  generalized  first  derivatives.  Assuming  that  po  on  Td  is  extended 
to  a  function  in  H1( fly),  multiplying  equation  (1)  by  a  test  function  q  £  Vf, 
equation  (3)  by  a  test  function  u  £  Ve,  and  integrating  by  parts,  we  obtain  the 
following  variational  form  of  (1)  -  (5):  Find  p  such  that  p  —  po  €  Vf,  and  u  G  Ve 
such  that  for  all  q  £  Vf  and  all  v  £  Ve, 


-Kf  +  k2Mf  -  ikGf 

-pfu2T 

P 

r 

rjp* 

— Ke  +  w2Me  _ 

u 

0 

—  J  VpS/q  +  k2  j  pq  —  ik  J  pq  —  J1  J  pf{v  ■  u)q  =  0, 

£if  f2  f  ra  r 

—  J  (A(V  •  u)(V  •  v)  +  2 pe(u)  :  e(v))  +  lu2  J  peu  •  v  —  J p(u  ■  v)  =  0. 

oe  r 

Proceeding  as  usual  in  the  Finite  Element  method  [2],  we  replace  Vf  and  Ve  with 
conforming  finite  element  spaces  and  obtain  the  algebraic  system 

(6) 

In  the  coupled  system  (6),  p  and  u  are  the  vectors  of  the  (values  of)  degrees  of 
freedom  of  p  and  u,  i.e.,  p  and  u  are  the  finite  element  interpolations  of  p  and  u, 
respectively.  The  matrix  blocks  in  (6)  are  defined  by 

P*K/q  =  /  vp-  Vq,  p*M/q  =  f  pq, 
nf  nf 

p*G/q  =  J  pq,  u*Kev  =  J  (A(V  •  u)(V  •  v)  +  2pe(u)  :  e(u)) , 
r„ 

u*Mev  =  Jpe(u-v),  p*Tv  =  J p(v  ■  v). 
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4.  SUBSTRUCTURED  DISCRETE  EQUATIONS 

The  fluid  and  solid  domains  are  decomposed  into  nonoverlapping  subdomains 
that  consist  of  union  of  elements, 


Nf 

%  =  UK 


Ne 


a 


=  u$- 


(7) 


The  vectors  of  degrees  of  freedom  corresponding  to  and  fl*  are  denoted  by  ps 
and  us,  respectively. 

The  subdomain  matrices  are  defined  by  subassembly, 


PS  K/q  =  /  Vp-Vg,  ps  M^qs  =  /  pq, 


ps*Gsfqs  =  I  pq , 
US*M*V  =  J  pe{u-v) 


J  pq ,  us*K*vs  =  J  A(V  •  u)(V  •  v)  +  2p,e(u)  :  e(v), 


pr*Trsvs  = 


J  p{w). 


We  will  use  vectors  consisting  of  all  subdomain  degrees  of  freedom, 


■  P1  ■ 

■  u1  ' 

p  = 

_pNf  _ 

,  u  = 

_u^  . 

and  the  corresponding  partitioned  matrices, 


. .  0 

'K  l  . 

..  O' 

K/  =  diag(K})  = 

,  Ke  =  diag(Kg)  = 

0  . 

l 

0  • 

KNe 

•  e  J 

The  matrices  M/ ,  Gf ,  and  Me  defined  similarly,  and 


T  =  (Trs) 


rpll  rj-il,7Ve 

rj^Nf,l  rj^Nf,Ne 


Let  N f  and  Ne  be  the  matrices  with  0, 1  entries  of  the  global  to  local  maps 
corresponding  to  the  decompositions  of  f If  and  fle,  respectively,  cf.,  (7),  so  that 

K/  =  N}K/N/,  Ke  =  N*KeNe. 

Let  B /  =  (Bj-, . . . ,  B^f )  and  Be  =  (Bj, . . . ,  B^e)  be  matrices  with  entries  0,  ±1 
such  that  the  conditions  B/p  =  0  and  Beu  =  0  express  the  constraint  that  the 
values  of  the  same  degrees  of  freedom  on  two  different  subdomains  coincide,  that 
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is. 


Byp  —  0  4=4>  p  =  Njp  for  some  p  (8) 

Beii  =  0  u  =  Neu  for  some  u.  (9) 

Multiplying  the  second  equation  in  (6)  by  to2pf  to  symmetrize  the  off-diagonal 
block  and  introducing  Lagrange  multipliers  A /  and  Ae  for  the  constraints  B/  p  =  0 
and  Beu  =  0,  we  get  the  system  of  linear  equations  in  block  form, 


'  -Kf  +  k2Mf  - 

ikG  —oj2pfT 

B  / 

0  ' 

■  p  " 

"  r  " 

LV2pf(-Ke+u2Me) 

0 

b: 

u 

0 

B  / 

0 

0 

0 

A/ 

0 

0 

Be 

0 

0 

.  Ae. 

.0. 

where  N*r  =  r.  The  system  (10)  is  equivalent  to  (6)  in  the  following  sense. 

Lemma  4.1.  If  (p,  u,  A/,  Ae)  is  a  solution  of  (10),  then  p  =  N/p  and  u  =  Neu, 
where  (p,u)  solves  (6).  Conversely,  z/(p,u)  is  a  solution  of  (6),  then  there  exist 
A /  and  \e  such  that  (p,u,  A/,  Ae)  is  a  solution  of  (10). 


Proof.  To  simplify  notation,  write  (6)  as  N*SNx  =  N*b  and  (10)  as 
Sx  +  B*A  =  b,  Bx  =  0.  Then  (8),  (9)  become  NullB  =  Range N. 

Suppose  that  Sx  +  B*A  =  b  and  Bx  =  0.  From  the  second  equation  and 
NullB  =  Range N,  it  follows  that  there  is  x  such  that  x  =  Nx.  Multiplying 
the  first  equation  by  N*  and  using  the  fact  that  N*B*  =  (BN)*  =  0,  we  obtain 
N*SNx  =  N*b. 

Conversely,  suppose  that  N*SNx  =  N*b.  Let  x  =  Nx.  Then  Bx  =  0,  and 
Sx  -  b  G  Null N*  =  (Range N)-1  =  (NullB)1-  =  Range B*, 


hence  Sx  —  b  =  B*A  for  some  A.  ■ 


Note  that  symmetry  of  the  system  or  the  existence  of  a  Lagrangean  functional 
are  not  needed. 

We  will  want  to  eventually  eliminate  the  variables  p  and  u.  Because  the  matrices 
— K/+fc2My  and  — Ke  +o>2Me  are  typically  close  to  singular  due  to  near-resonance, 
we  replace  them  by  regularized  matrices 

A  j  =  —  Ky  +  fc2My  +  ikGf  +  R/, 

Ae  =  — Ke  +  W2  Me  +  R-e, 


where  the  regularization  matrices  are  given  by 


R/  =  diag(R')), 
Re  =  diag(R*), 


N  , 


ps*R)qfi  =  ik^, 


m. 


t= l  ° 

t^s 

Ne 


us*R*vs  =  iujpe  i 


( n  ■  u)(n  ■  v), 


t^s  dQ.%r\dO.% 
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where  of  G  {0,  ±1},  of  =  - of ,  and  of  =  {0,  ±1},  of  =  -a*/. 

It  is  shown  in  [6]  that  if  for  a  given  s,  all  of  >  0  or  all  of  <  0  with  some 
of  ^  0,  then  A).  is  invertible.  The  case  of  solid  subdomains  is  similar.  For  details 
on  strategies  for  choosing  of  to  guarantee  this,  see  [6].  In  our  computations,  we 
simply  choose  of  =  +1  if  s  >  t,  of  =  —  1  if  s  <  t,  and  similarly  for  of 
Because 

N}R/N/=0,  N*ReNe  =  0, 
the  system  (10)  is  equivalent  to 


A/ 

-w2p/T 

B/ 

0  ' 

"  p  ' 

"  r  " 

-u2PfT* 

u)2pfAe 

0 

b: 

u 

0 

B  / 

0 

0 

0 

A/ 

0 

0 

Be 

0 

0 

.  Ae. 

.  0. 

Since  the  value  of  Tu  depends  on  the  values  of  u  on  T  only,  we  have 

Tu  =  TJeur,  ur  =  JgU, 

where  Je  is  the  matrix  of  the  operator  of  embedding  a  subvector  that  corresponds 
to  degrees  of  freedom  on  T  into  u  by  adding  zero  entries.  Similarly, 

T*p  =  T*Jypr,  pr=J}p. 

Therefore,  we  obtain  the  augmented  system  equivalent  to  (11), 


"A  / 

0 

B  / 

0 

0 

~co2pfTJe 

p 

r 

0 

u2pfAe 

0 

b: 

—u>2pfT*Jf 

0 

u 

0 

B  / 

0 

0 

0 

0 

0 

A  / 

0 

0 

Be 

0 

0 

0 

0 

Ae 

0 

T* 

J/ 

0 

0 

0 

-I 

0 

Pr 

0 

0 

T* 

Je 

0 

0 

0 

-I 

ur 

0 

Because  the  variables  in  a  coupled  system  typically  have  vastly  different  scales, 
we  use  symmetric  diagonal  scaling  to  get  the  scaled  system 


'A/ 

0 

B  / 

0 

0 

-TJe' 

P 

r 

0 

Ae 

0 

b: 

— T*J/ 

0 

u 

0 

B  / 

0 

0 

0 

0 

0 

A  / 

0 

0 

Be 

0 

0 

0 

0 

Ae 

0 

T* 

J/ 

0 

0 

0 

-I 

0 

Pr 

0 

0 

T* 

Je 

0 

0 

0 

I 

ur 

0 

where  the  matrices  and  the  vectors  scale  as 

Af  =  D/A/D/,  Ae  =  u!2  pfDeAeDe,  T  =  w2p/D/TDe,  (14) 
B  /  =  E/B/D/,  Be  =  EeBeDe,  r  =  D/f,  (15) 

p  =  D/p,  u  =  Deu,  A/ =  D/A/,  Ae  =  DeAe.  (16) 
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The  scaling  matrices  Dy,  De,  Ey,  and  Ee,  are  diagonal.  We  have  chosen  scaling 
matrices  with  positive  diagonal  entries  such  that  the  absolute  values  of  the  diagonal 
entries  of  A y  and  Ae  are  one  and  the  £2  norms  of  the  columns  of  Be  and  B y  are 
one.  Other,  perhaps  better,  scaling  choices  are  of  course  possible. 

The  proposed  method  consists  of  eliminating  the  variables  p  and  u  from  the 
augmented  system  (13)  and  solving  the  resulting  reduced  system  by  Generalized 
Conjugate  Residuals,  with  preconditioning  by  projection  on  a  coarse  space  with  few 
degrees  of  freedom  per  subdomain. 

Computing  p  and  u  from  the  first  two  equations  in  (13)  gives 

p  =  A71(f-B}A/  +  TJeur)  (17) 

u  =  A“1(— B*Ae  +  T*J/pr)  (18) 


Substituting  p  and  u  from  (17),  (18)  into  the  rest  of  the  equations  in  (13),  we 
obtain  the  reduced  system 

Fx  =  b,  (19) 

where 


ByA/B} 

0 

0  -ByA^TJ 

0 

BeA- 

!Be  — BeAj1T*Jy 

T*  A  _1R* 

JfAf  Uy 

0 

-I  JJA^TJe 

0 

— Je  Ae 

4Be  JeA”lrr*Jy  -I 

and 


i 

lU 

1  ^ 

'< 

i _ 

Ae 

U  _ 

0 

X  — 

Pr 

,  D  — 

.  ur  . 

0 

In  equation  (20),  the  first  diagonal  block  ByAj^By  is  exactly  same  as  in 
the  FETI-H  method  for  the  Helmholtz  equation.  The  second  diagonal  block 
ByA^B}  is  the  analogue  of  FETI-H  for  the  elastodynamic  problem.  FETI-H 
is  known  to  converge  fast,  so  one  can  expect  that  on  the  subspace  defined  by 
the  coarse  correction  (cf.,  Section  5  below),  these  two  diagonal  blocks  will  be  well 
conditioned.  The  off-diagonal  blocks  all  contain  the  inverse  of  Ay,  resp.  Ae, 
which  are  discretizations  of  differential  operators.  Hence,  the  off-diagonal  blocks 
are  discretizations  of  compact  operators,  and  the  structure  of  the  equation  (19) 
resembles  discretization  of  a  Fredholm  integral  equation  of  the  second  kind,  which 
is  well  conditioned  because  its  eigenvalues  cluster  about  a  point  different  from  zero. 
Such  clustering  is  indeed  observed  computationally,  cf.,  Section  7  below. 

Evaluating  the  matrix  vector  product  Fx  requires  the  solution  of  one  independent 
problem  per  subdomain,  because 


’  xf  ’ 

Byq 

Ae 

Bev 

f 

F 

Pr 

— 

Jy9l  -  Pr 

,  where  < 

.  ur  . 

.  J*v  -  ur  . 

q  =  Ay  1(  — ByAy  +  TJeur), 
v  =  A“!(— B*Ae  +  T*Jypr). 
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5.  ITERATIVE  SOLUTION 

The  present  method  consists  of  solving  the  linear  system  (19)  by  the  Generalized 
Conjugate  Residual  (GCR)  method  [12]  preconditioned  by  a  subspace  correction. 
GCR  is  mathematically  equivalent  to  GMRES  and  easier  to  program,  but  less 
numerically  stable.  Because  the  number  of  iterations  was  small,  GCR  was  sufficient 
for  our  purposes. 

Let  Q  be  a  matrix  with  the  same  number  of  rows  as  F  and  linearly  independent 
columns.  The  columns  of  Q  form  the  basis  of  the  coarse  space.  We  wish  to  enforce 
the  condition  that  the  residual  is  orthogonal  to  the  coarse  space, 


Q*(Fx  —  b)  =  0,  (21) 

by  adding  a  correction  from  the  coarse  space  in  each  iteration.  Define  the  correction 
operator  C  by 


C(v,b)  =  v  +  Qw,  Q*(F(v  +  Qw)  -  b)  =  0. 

We  first  make  make  sure  that  the  initial  approximation  satisfies  (21):  If  v  is  given, 
the  initial  approximation  is  =  C(v,b).  The  system  (19)  is  then  solved  by  the 
GCR  method  with  left  preconditioning  by  the  action  of  the  projection  P,  given  by 

Pv  =  C(v,  0)  =  (I  -  Q(Q*FQ)_1Q*F)v. 

Equivalently,  the  GCR  method  is  applied  to  the  preconditioned  system 


PFx  =  Pb.  (22) 

In  step  to,  the  GCR  method  computes  x("0  by  adding  to  a  linear 

combination  of  preconditioned  residuals  P(Fx^fc)  —  b),  k  =  0, . . . ,  to  —  1,  so  that  the 
£2  norm  of  the  residual  ||P(Fx("0  —  b)||  is  minimal.  Because  Q*FP  =  0,  it  follows 
that  the  iterates  run  in  a  subspace, 

Q*F(x(m)  _x(0))  =0, 

hence,  from  the  selection  of  x^°\ 

Q*(Fx(m)  -  b)  =  0. 

We  choose  the  matrix  Q  of  the  form 

'  D/B/diag(Y*)s  0  0  0 

Q  =  0  DeBediag(Y®)s  0  0 

^  0  0  D/J)diag(Z})s  0 

0  0  0  DeJ*diag(Z*)s . 

For  a  fluid  subdomain  12  J,  we  choose  Y)  as  the  matrix  of  columns  that  are  discrete 
representations  of  plane  waves  in  a  small  number  of  equally  distributed  directions, 
or  discrete  representation  of  the  constant  function.  For  a  solid  subdomain  12®,  the 
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FIG.  1.  Model  2D  Problem 

r. 


columns  of  Y|  are  discrete  representations  of  plane  pressure  and  shear  waves,  or 
of  rigid  body  motions.  The  matrices  Zsj  and  Zse  are  chosen  in  the  same  way  as  Y  j 
and  Y|,  with  possibly  different  selection  of  the  number  of  directions  and  selection 
of  constant  or  rigid  body  modes.  Some  of  the  matrices  Yj,  Y|,  Z®,  or  Y|  may  be 
void. 


6.  COMPUTATIONAL  RESULTS 

We  consider  a  model  2D  problem  with  a  scatterer  in  the  center  a  waveguide,  cf., 
Fig.  1.  The  fluid  domain  flf  is  a  square  with  side  1  in,  filled  with  water  with  density 
Pf  =  1000  kgm~3  and  speed  of  sound  c/  =  1500  ?ns-1.  The  scatterer  is  a  square 
in  the  center  o  the  fluid  domain,  with  side  0.2  m,  and  consisting  of  aluminum  with 
dnsity  pe  =  2700  kg  m~3  and  Lame  elasticity  coefficients  A  =  5.5263.10 10  Nm~2, 
p.  =  2.595.10 10Nm~2. 

For  numerical  stability,  the  rows  of  the  matrices  Be,  B/,  and  the  columns  of  Q 
were  orthogonalized.  The  number  of  reduced  variables  reported  in  the  tables  is  the 
number  of  the  Lagrange  multipliers  to  enforce  the  linearly  independent  constraints. 
The  number  of  coarse  degrees  of  freedom  is  the  size  of  the  basis  of  the  columns  of 
the  matrix  Q. 

The  fluid  domain  and  the  solid  domain  are  discretized  by  bilinear  square  elements 
on  a  uniform  mesh  with  meshsize  h.  Both  domains  are  divided  independently  into 
m  by  n  subdomains  by  dividing  their  horizontal  sides  into  m  intervals  of  the  same 
length  and  the  vertical  sides  into  n  intervals  of  the  same  length.  The  iterations 
were  stopped  when  the  relative  residual  in  the  £2  norm  of  the  reduced  and  scaled 
system  (19)  was  less  than  10-6. 

The  numerical  results  for  the  model  problem  are  summarized  in  Tables  1  to  3. 
The  lines  in  the  tables  are  organized  in  groups.  The  first  line  in  the  group  is 
the  result  for  the  FETI-H  method  with  Neumann  condition  instead  of  the  elastic 
scatterer.  The  following  lines  are  for  the  elastic  scatterer,  and  they  show  the  number 
of  iterations  for  increasing  the  number  of  coarse  space  functions,  that  is,  the  columns 
of  Q.  In  all  cases,  the  constant  function  is  included  in  the  coarse  space  for  the  fluid 
and  the  rigid  body  modes  are  included  in  the  coarse  space  for  the  solid,  both  for 
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TABLE  1 


Decreasing  mesh  size  and  wavelength, 


Problem  description 

Coarse  directions 

Subdomains 

Multipliers 

1  Wet  interf. 

h 

k 

Fluid 

Solid 

FI. 

Sol. 

FI. 

Sol. 

1/40 

10 

4x4 

rigid 

4 

1/40 

10 

4x4 

lxl 

4 

0 

0 

0 

1/40 

10 

4x4 

lxl 

4 

0 

4 

0 

1  solid  subdomain 

Number  of 
degrees  of  freedom 


Orig. 

Red. 

Coarse 

1632 

240 

63 

1794 

312 

63 

1794 

312 

86 

Iter. 


11 

18 

13 


1/80 

20 

4x4 

rigid 

4 

6336 

464 

63 

1/80 

20 

4x4 

lxl 

4 

0 

0 

0 

6914 

600 

63 

1/80 

20 

4x4 

lxl 

4 

0 

4 

0 

6914 

600 

86 

1/160 

40 

4x4 

rigid 

4 

24960 

912 

63 

1/160 

40 

4x4 

lxl 

4 

0 

0 

0 

27138 

1176 

63 

1/160 

40 

4x4 

lxl 

4 

0 

4 

0 

27138 

1176 

86 

1/160 

40 

4x4 

rigid 

8 

24960 

912 

102 

1/160 

40 

4x4 

lxl 

8 

0 

4 

0 

27138 

1176 

125 

1/160 

40 

4x4 

lxl 

8 

4 

4 

4 

27138 

1176 

126 

TABLE  2 


Increasing  number  of  solid  subdomains 


Problem  description 

Coarse  directions 

Number  of 

Subdomains 

Multipliers 

Wet  interf. 

degrees  of  freedom 

h 

k 

Fluid 

Solid 

FI. 

Sol. 

FI. 

Sol. 

Orig. 

Red. 

Coarse 

Iter. 

1/250 

20 

4x4 

rigid 

4 

60600 

1416 

63 

32 

1/250 

20 

4x4 

lxl 

4 

0 

0 

0 

65802 

1824 

63 

37 

1/250 

20 

4x4 

lxl 

4 

0 

4 

0 

65802 

1824 

86 

32 

1/250 

20 

4x4 

rigid 

4 

60600 

1416 

63 

32 

1/250 

20 

4x4 

2x2 

4 

0 

0 

0 

65802 

2034 

72 

48 

1/250 

20 

4x4 

2x2 

4 

0 

4 

0 

65802 

2034 

100 

42 

1/250 

20 

4x4 

rigid 

4 

60600 

1416 

63 

32 

1/250 

20 

4x4 

3x3 

4 

0 

0 

0 

65802 

1824 

63 

37 

1/250 

20 

4x4 

3x3 

4 

0 

4 

0 

65802 

1824 

86 

32 

1/250 

20 

4x4 

rigid 

8 

60600 

1416 

102 

26 

1/250 

20 

4x4 

lxl 

8 

0 

4 

0 

65802 

1824 

125 

27 

1/250 

20 

4x4 

lxl 

8 

4 

4 

4 

65802 

1824 

126 

27 

1/250 

20 

4x4 

rigid 

8 

60600 

1416 

102 

26 

1/250 

20 

4x4 

2x2 

8 

0 

4 

0 

65802 

2034 

139 

39 

1/250 

20 

4x4 

2x2 

8 

4 

4 

4 

65802 

2034 

148 

36 

1/250 

20 

4x4 

rigid 

8 

60600 

1416 

102 

26 

1/250 

20 

4x4 

3x3 

8 

0 

4 

0 

65802 

1824 

125 

27 

1/250 

20 

4x4 

3x3 

8 

4 

4 

4 

65802 

1824 

126 

27 
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TABLE  3 

Increasing  number  of  fluid  and  solid  subdomains 


Problem  description 

I  Subdomains 


Coarse  directions 
Multipliers  I  Wet  interf. 


Number  of 
degrees  of  freedom 


the  multipliers  and  for  the  wet  interface.  The  prototype  implementation  was  done 
in  MATLAB,  hence  we  do  not  report  timings. 

Table  1  shows  that  the  method  is  scalable  when  the  frequency  is  increased  and  the 
mesh  is  correspondingly  refined.  Since  there  is  only  one  solid  subdomain,  the  effect 
of  increasing  the  coarse  space  for  the  solid  and  for  the  wet  interface  was  minor.  The 
results  in  Table  2  show  that  the  number  of  iterations  does  not  increase  when  the 
scattered  is  decomposed  into  increasing  number  of  subdomains.  Finally,  Table  3 
shows  the  number  of  iterations  does  not  grow  when  a  fixed  problem  is  decomposed 
into  a  growing  number  of  both  fluid  and  solid  subdomains. 

7.  SPECTRUM  OF  THE  ITERATION  OPERATOR 

Convergence  of  the  method  can  be  explained  by  an  analysis  of  the  spectrum  of 
the  iteration  operator.  From  well-known  properties  of  Krylov  space  methods  [12],  it 
follows  that  in  m  iterations,  the  residual  is  reduced  by  the  factor  of  at  least  Cr(m), 
where 

r(m)  =  min  max  |»(A)|. 

deg(p)=TO  Jvep(PF) 
p(0)  =  l 

Here,  the  minimum  is  taken  over  complex  polynomials  and  the  maximum  over 
eigenvalues  of  PF. 

Figures  2  and  3  show  one  representative  case  of  the  spectrum  of  the  iteration 
operator  PF  from  (22)  for  Neumann  boundary  condition  on  the  scatterer,  and  for 
the  elastic  scatterer,  respectively.  The  problem  setting  was  h  =  1/40,  A:  =  10,  4  x  4 
fluid  subdomains  and,  for  Figure  3,  2  x  2  solid  subdomains.  The  coarse  space  had 
4  directions  for  fluid  waves,  both  for  the  multipliers  and  for  the  wet  interface. 

The  figures  show  that  in  both  cases,  the  spectrum  is  clustered  around  a  point 
other  than  the  origin.  For  the  coupled  problem,  there  are  few  more  eigenvalues  near 
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FIG.  2.  Spectrum  of  the  FETI-H  iteration  operator 


Complete  spectrum  Spectrum  detail 


FIG.  3.  Spectrum  of  the  iteration  operator  for  coupled  problem 

Complete  spectrum  Spectrum  detail 


the  origin  and  some  clustered  around  —1.  Because  the  eigenvalues  are  clustered 
around  a  point  different  from  zero,  a  polynomial  with  a  low  value  of  r(?n)  exists, 
and  convergence  is  fast.  The  effect  of  adding  the  elastic  scatterer  shows  in  few  extra 
eigenvalues  close  to  the  cluster  and  a  minor  change  in  other  eigenvalues. 

The  change  in  the  spectrum  can  be  heuristically  explained  as  follows,  cf.,  also 
the  remarks  at  the  end  of  Section  5.  The  last  two  by  two  diagonal  block  (with 
—I  on  the  diagonal)  in  (20)  is  a  discretization  of  a  Fredholm  integral  equation  of 
the  second  kind,  with  eigenvalues  clustered  around  —  1.  The  influence  of  the  other 
off-diagonal  blocks  is  relatively  weak.  The  FETI-H  method  already  has  eigenvalues 
around  —1,  so  the  overall  effect  of  adding  the  elastic  scatterer  on  the  spectrum  of 
the  iteration  operator  is  small. 

8.  CONCLUSION 

We  have  presented  a  new  fast  substructuring  method  for  coupled  fluid-solid 
acoustic  problems.  The  computational  results  indicate  that  the  method  is  scalable 
with  respect  to  mesh  size,  frequency,  and  the  number  of  subdomains.  The  growth 
of  the  number  of  iterations  can  be  controlled  by  increasing  the  size  of  the  coarse 
space.  Numerical  calculation  of  the  spectrum  of  the  iteration  operator  suggests 
that  the  fast  convergence  is  due  to  clustering  of  the  spectrum.  In  most  cases,  the 
resulting  number  of  iterations  is  same  or  slightly  larger  than  for  the  FETI-H  method 
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for  the  same  problem  with  the  Neumann  boundary  condition  instead  of  an  elastic 
scatterer. 

Performance  study  of  the  method  for  realistic  problems,  3D  computations, 
parallel  implementation,  and  theoretical  analysis  will  be  reported  elsewhere. 
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